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ABSTRACT. In the general context of complex data processing, this paper 
reviews a recent practical approach to the continuous wavelet formalism on 
the sphere. This formalism notably yields a correspondence principle which 
relates wavelets on the plane and on the sphere. Two fast algorithms are also 
presented for the analysis of signals on the sphere with steerable wavelets. 



1. Introduction 

In many fields of science, from computer vision, to biomedical imaging, 
geophysics, or astrophysics and cosmology, experiments are set up releasing 
more and more complex data to process. A first complexity of the data 
lies in their large volume, related to the always increasing resolution of 
technological devices. Moreover, data are not necessarily distributed on the 
real line (audio signals, ...), or on the plane (images, ...), but can live on 
higher-dimensional or nontrivial manifolds (W 1 , sphere, hyperboloid, ...). 
Finally, the data may correspond not only to scalar fields (local intensity), 
but also to tensor fields on those manifolds (local diffusion matrix, local 
polarization, ...). 

In this new era of complex data processing, powerful tools always need 
to be developed for the precise analysis of the signals under scrutiny. In this 
paper, we review recent formal and algorithmic advances for the continuous 
wavelet analysis of signals on the sphere. This scale-space formalism goes 
well beyond the spectral analysis, as it enables one to probe the localization, 
scale, and orientation of the features of the signals analyzed. An exhaus- 
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tive review on wavelets on the sphere and related manifolds is presented in 
another article of the present issue [5] . 

For the sake of the illustration we choose the example of the cosmic 
microwave background (CMB) data from cosmology. The CMB is a polar- 
ized electromagnetic blackbody radiation observed today in all directions 
of the sky, which emerged some 380.000 years after the Big Bang. This 
snapshot of the early universe bears a wealth of information for the study of 
its structure and evolution, i.e., for cosmology. The present NASA WMAP 
( Wilkinson Microwave Anisotropy Probe) satellite experiment [6] releases 
maps of the celestial sphere of 3 megapixels at each detection frequency, 
while the forthcoming ESA Planck Surveyor satellite experiment |24j will 
increase the resolution to 50 megapixels. The CMB therefore crystallises 
the previously quoted potential data complexities. Its temperature (inten- 
sity) and polarization, respectively, define scalar and tensor fields on the 
sphere, and the corresponding experimental data already appear in large 
volumes. Various applications of the continuous wavelet formalism on the 
sphere for the analysis of the CMB are presented in another article of the 
present issue [23] . 

The structure of the paper goes as follows. We only focus on the for- 
malism for the continuous wavelet transform on the sphere introduced in [2 J , 
as recently further developed in a practical approach by [29]. This formalism 
is explicitly reviewed in section [2j The wavelet decomposition of a signal on 
the sphere S 2 is defined by its projection coefficients on translated, rotated, 
and dilated versions of a mother wavelet, i.e., a directional correlation at 
each analysis scale. These wavelet coefficients therefore live on the rotation 
group in three dimensions 50(3). The wavelet must satisfy an admissibility 
condition ensuring that the signal may be explicitly reconstructed from its 
wavelet coefficients. A correspondence principle is also recalled stating that 
wavelets on the sphere may be built from an inverse stereographic projec- 
tion of wavelets on the plane. This principle enables one to transfer onto 
the sphere some properties of wavelets on the plane, such as the notion of 
steerability. We explicitly describe major examples of axisymmetric, direc- 
tional, and steer able wavelets. In section El we give the generic definition of 
directional correlation, and the definition of standard correlation, to which 
reduce the wavelet coefficients of a signal with a steerable or axisymmetric 
wavelet. We discuss their a priori computation cost on any pixelization of 
S 2 and of 50(3), which is prohibitive for high resolution data. We em- 
phasize the existence of a directional and standard correlation relation in 
harmonic space. In section U we first discuss the band- limitation of signals 
and filters. We then review two fast algorithms for the directional correla- 
tion of band-limited signals and filters on iso-latitude pixelizations on the 
sphere. The first one is based on a technique of separation of variables in 
the Wigner .D-functions on 50(3) [18] [30]. The second one relies on the fac- 
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torization of the three-dimensional rotation operators to interpret the result 
of the directional correlation as a function on the three-torus T 3 , and ap- 
plies the separation of variables to three-dimensional imaginary exponentials 
\25 \ \28 \ [22] . The a priori 0(L 5 ) asymptotic complexity is thereby reduced 
to 0(L 4 ), where 2L roughly stands for the square-root of the number of 
pixels on the sphere, i.e., for band- limited signals and filters with band-limit 
L £ N. For steerable and axisymmetric wavelets, the directional correla- 
tion resumes to standard correlations, and the asymptotic complexity drops 
to 0(L 3 ). The typical computation time for the directional correlation of 
megapixels maps (L ~ 10 3 ) correspondingly drops from years to tens of 
seconds on a standard computer. This easily allows the analysis of multiple 
signals at such high resolutions, and at multiple scales. These developments 
finally lead us to our conclusions in section [5j 

2. Continuous wavelets on the sphere 
2.1 Practical approach 

Among other approaches |17} [TT1 [12] , a satisfactory formalism for the contin- 
uous wavelet transform of signals on the sphere S 2 was originally established 
in a group-theoretical framework [21 EJ HJ El Ej . The aim of the present article 
is to review a more practical but completely equivalent approach, recently 
proposed by [29]. In that framework, a "mother wavelet" ty(uj) is defined as 
a localized square-integrable function on the unit sphere, on which contin- 
uous affine transformations (translations, rotations, and dilations) may be 
applied. The wavelet transform of a square-integrable signal on the sphere 
is then defined as the directional correlation of the signal with the dilated 
versions of the mother wavelet. At each scale, the corresponding wavelet 
coefficients are square-integrable functions on the rotation group in three 
dimensions SO (3). Finally, an admissibility condition is imposed on the 
wavelet which ensures an exact reconstruction formula of the signal from its 
wavelet coefficients . 

The real and harmonic structures of the unit sphere S 2 are concisely 



Notice that the signals and filters considered by the formalism are scalar functions, 
i.e., invariant under local rotations in the tangent plane at each point on the sphere. 
In the general context of complex data processing, one might want to generalize 
the wavelet formalism presented here to the analysis of rank n tensor functions. 
However, tensor fields may equivalently be expressed in terms of scalar fields. In 
particular, polarization data on the sphere constitute a rank 2 tensor field. It 
can be equivalently described in terms of its so-called electric and magnetic parts, 
which actually constitute two separate functions on the sphere, with a purely scalar 
behaviour. The present formalism for the scale-space wavelet decomposition of 
scalar functions on the sphere may therefore be applied to the analysis of both 
scalar fields, and tensor fields such as polarization data |30l 131] . 



4 



Y. Wiaux, J. D. McEwen, and P. Viclva 



summarized as follows. Any point uj on the sphere is given in spherical 
coordinates asw = (9, ip), in terms of a polar angle, or co-latitude 6 £ [0, tt], 
and an azimuthal, or longitudinal angle ip £ [0,2w[. Let G(u) be a square- 
integrable function on the sphere, i.e., G(uS) in L 2 (S 2 , dd), with the invariant 
measure d£l = d cos 6 dip. The spherical harmonics form an orthonormal 
basis for the decomposition of functions in L 2 (S 2 ,dQ,). They are explicitly 
given in a factorized form in terms of the associated Legendre polynomials 
P[ n (cos9) and the complex exponentials e imif as 



Y lm {6,<p) 



2l + l(l-m)\ 
47r (l + m)\ 



1/2 

P{ n (cos9)e im ^, (2.1) 



with I £ N, m E Z, and |m| < / [TJ [27] . While the index Z represents 
an overall frequency on the sphere, \m\ represents the frequency associated 
with the azimuthal variable (p. Any G(uS) is thus uniquely given as a linear 
combination of scalar spherical harmonics G (uj) = X)zgn S| m |<z Gi m Yi m (uj) 
(inverse transform), for the scalar spherical harmonics coefficients G[ m = 
j S 2 d^lY^ (uj) G (uj) (direct transform), with \m\ < I. 

The continuous affine transformations on functions on the sphere are 
defined as follows. The operator R(ujq) for the motion, or translation, of 
amplitude ujq = (9q, ipo) of a function reads 

[R(uj )G](uj) = G{R^uj), (2.2) 

where R Uo (9,(p) = [R^ q Rq ](9, ip) is defined by the three-dimensional rota- 
tion matrices Rg and Rp , acting on the Cartesian coordinates (x,y,z) in 
three dimensions centered on the sphere and associated with uj = (9, cp). The 
rotation operator R z (x) of a function around itself, by an angle x G [0) 27r[, 
is given as 

R s ( X )g] (uj) = g(r^ V), (2.3) 

where R^(9,ip) = (9,<p + x) a ls° follows from the action of the three- 
dimensional rotation matrix R* on the Cartesian coordinates (x, y, z) associ- 
ated with uj = (9, <p). The dilation operator D(a) on functions in L 2 (S 2 , dQ), 
for a dilation factor a £ , is defined in terms of the inverse of the corre- 
sponding dilation D a on points in S 2 as 

[D (a) G] (uj) = A 1 / 2 (a, 9) G (L>~M , (2-4) 

with X 1 / 2 (a,9) = a _1 [l + tan 2 (#/2)]/[l + a' 2 tan 2 (6>/2)]. The dilated point 
is given by D a (9,ip) = (9 a (9),(p) with the linear relation tan(0 a (0)/2) = 
atan(0/2). The dilation operator therefore maps the sphere without its 
South pole on itself: 9 a (9) : 9 £ [0,ir[— ► 9 a £ [0, w[. This dilation operator is 
uniquely defined by the requirement of the following natural properties |29] . 
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The dilation D a of points on S 2 must be a radial (i.e., only affecting the 
radial variable 9 independently of tp, and leaving p invariant) and conformal 
(i.e., preserving the measure of angles in the tangent plane at each point 
of S 2 ) diffeomorphism (i.e., a continuously differentiable bijection on S 2 ). 
The factor A(a, 9) explicitly appears in the conformal transformation of the 
metric through the dilation D a . The normalization by X 1 ^ 2 (a,9) in (|2.4h is 
uniquely determined by the requirement that the dilation D{a) of functions 
in L 2 (S 2 ,dQ) be a unitary operator (i.e., preserving the scalar product in 
L 2 (S 2 ,dQ), and specifically the norm of functions). 

The analysis of signals goes as follows. The wavelet transform of a 
signal F(uj) in L 2 (S 2 ,dQ) on the sphere, with the wavelet ^(co), localized 
analysis function in L 2 (S 2 ,d£l), is defined as the directional correlation be- 
tween F(u) and the dilated wavelet "if a = D(a)^, i.e., as the scalar product: 

w£ ( P , a) = (* p>a \F) = I dn % >a h f m , (2.5) 

J 

with ^ Pt a = R(p)^a, and p = (Oo,cpo,x)- At each scale, the wavelet coef- 
ficients W£(p, a) are therefore square-integrable functions on the rotation 
group in three dimensions SO (3). They represent the characteristics of the 
signal for each analysis scale a, direction x-> an d position uq. This defines 
the scale-space nature of the wavelet decomposition on the sphere. 

The real and harmonic structures of the rotation group in three dimen- 
sions 50(3) are concisely summarized as follows. Any rotation p on S0(3) 
is given in terms of the three Euler angles p = (9,p,x), with 9 £ [0, tt], and 
ip, x £ [0, 27r[. Let H (p) be a square-integrable function on 50(3), i.e., H(p) 
in L 2 (50(3), dp), with the invariant measure dp = dipd cos 9dx- The Wigner 
D-functions are the matrix elements of the irreducible unitary representa- 
tions of weight I of the group in L 2 (SO(3),dp). By the Peter- Weyl theorem 
on compact groups, the matrix elements also form an orthogonal basis 
in L 2 (SO(3),dp). They are explicitly given in a factorized form in terms 
of the real Wigner d-functions d l mn {9) and the complex exponentials, e~ tmLp 
and 6 _ as 

D l mn 9, x) = e- im *d l mn (9) e"** (2.6) 

with leN, m,fieZ, and |m|, \n\ < I |271 [8]. Again, I represents an overall 
frequency on 50(3), and \m\ and \n\ the frequencies associated with the 
variables p and x, respectively. Any H(p), such as the wavelet coefficients 
at each scale of a signal on S 2 , is thus uniquely given as a linear combination 
of Wigner ^-functions as H (p) = £ ZeN (2Z + 1)/8tt 2 E|m|,|n|<J^mn^mn (p) 
(inverse transform), with, for |m|,|n| < I, the Wigner -D-functions coeffi- 
cients H l mn = f SO / 3 \ dpD l mn (p) H (p) (direct transform). 

The synthesis of a signal F(u) from its wavelet coefficients reads as: 
r+oo An r 

F(u)= ^ dpW£(p,a)[R(p)L^ a ](Lu). (2.7) 

Jo ^ J so{3) 
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In this relation, the operator in L 2 (S 2 ,dQ) is defined □ by the follow- 
ing action on the spherical harmonics coefficients of functions: L^Gi m = 
Gi m /C\,, with |m| < I. This exact reconstruction formula holds if and only 
if the spherical harmonics transform \&/ m of the wavelet *$>{uS) satisfies the 
following admissibility condition [29J: 




for all I £ N. This condition intuitively requires that the whole wavelet 
family ^f a (cj), for a £ M5_, covers each frequency index I with a finite and 
non-zero amplitude. As explicitly expressed in section [31 the direct Wigner 
D-functions transform of the wavelet coefficients of a signal F with ^ is given 
as the pointwise product of the spherical harmonics coefficients i*) m and 



(^a)i n - The admissibility condition consequently requires that the wavelet 
family as a whole preserves the signal information at each frequency I. 

2.2 Correspondence principle 

Wavelets on the plane are well-known, and may be easily constructed as 
the corresponding admissibility condition reduces to a zero-mean condition 
for a function both integrable and square-integrable. On the contrary, the 
admissibility condition (|2.8p for wavelets on the sphere is difficult to check 
in practice. In that context, a correspondence principle was proved in |29j . 
stating that the inverse stereographic projection of a wavelet on the plane 
leads to a wavelet on the sphere. 

The stereographic projection is the unique radial conformal diffeomor- 
phism mapping the sphere S 2 onto the plane M 2 . The unitary stereographic 
projection operator between functions G in L 2 (S 2 ,dfl) and g in L 2 (R 2 , d 2 x), 
and its inverse, respectively read 



in spherical coordinates on the sphere u = (9,<p), and polar coordinates on 
the plane x = (r, </>) . The azimuthal coordinates on the plane and on the 
sphere are identified to one another: ip. The radial conformal diffeomor- 
phism between points is given as 7r(6,(p) = (r(9),ip) for r(6) = 2tan(6*/2), 



2 The operator in our notations coincides with the inverse of the standard frame 
operator A* defined in [5]- 





(2.8) 




(2.9) 
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and its inverse reads 7r _1 (r, ip) = (9(r),ip) for 9(r) = 2 arctan(r/2). The 
diffeomorphism r{6) and its inverse 9{r) explicitly define the stereographic 
projection and its inverse. This stereographic projection maps the sphere, 
without its South pole, on the entire plane: r(9) : 9 £ [0,ir[-^ [0, oof. Geo- 
metrically, it projects a point u> = (9, ip) on the sphere onto a point x = (r, ip) 
on the tangent plane at the North pole, co-linear with u and the South pole 
(see Fig. [1]). The pre-factors in (|2,9p are required to ensure the unitarity of 
the projection operators II and U . 

In this framework, the correspondence principle established states that, 
if the function Tp(r, p) in L 2 (M 2 , d 2 x) satisfies the wavelet admissibility con- 
dition on the plane, i.e., essentially a zero- mean condition, then the function 



in L 2 (S 2 , d£l), satisfies the wavelet admissibility condition (|2.8|) on the sphere. 
This enables the construction of wavelets on the sphere by projection of 
wavelets on the plane. It also transfers wavelet properties from the plane 
onto the sphere, such as the steerability discussed in the next subsection. 



FIGURE 1: Stereographic projection it and its inverse tt~ , relating points 
(9, ip) on the sphere and (r, p) on its tangent plane at the North pole. The 
same relation holds through IT and IT -1 between functions living on each 
of the two manifolds, as illustrated by the shadow on the sphere and the 
localized region on the plane |29j. 



2.3 Axisymmetric, directional, and steerable wavelets 

We present here axisymmetric, directional, and steerable wavelets on the 
sphere, built as inverse stereographic projections of wavelets on the plane. 



*(0,¥>) = [n-ty] (e,<p) 



(2.10) 




s 
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An axisymmetric filter is by definition invariant under rotation around 
itself. That is, when located at the North pole, an axisymmetric filter is de- 
fined by a function A{6) independent of the azimuthal angle <p. On the plane, 
the Mexican hat wavelet is defined as the normalized (negative) Laplacian 
of a Gaussian e^ x +v >f 2 . Its inverse stereographic projection defines the 
Mexican hat wavelet on the sphere (see Fig. [2]). 

Any non-axisymmetric filter is said to be directional, and is given as a 
general function ty(9,(p) in L 2 (S 2 ,dQ,). The elliptical Mexican hat wavelet 
is a directional modification of the axisymmetric Mexican hat, obtained by 
considering different widths a x and a y , respectively in the x and y directions 
on the plane for the original Gaussian: e^+^Z 2 -► e (* 2 /^+?/ 2 /^)/2 ^ 
The wavelet obtained as inverse stereographic projection of the normalized 
(negative) Laplacian of this Gaussian reads (see Fig. [2]) as: 



^ (mex) H = ^N(a x ,a y ) (l + tan 2 £) 



4 tan 2 (9/2 

CT 2 +CT 2 



cos 2 ip H — §■ sin 2 <p 



e 



-2tan 2 | (cos 2 ^/cr 2 +sin 2 </V°y) 



(2.11) 



The constant N(a x , a y ) = (a 2 + a 2 ) [a x a y (3a x + 3cr^ + 2ct 2 ct 2 )/2]~ 1//2 stands 
for the normalization. One can identify the wavelet parameters through the 
eccentricity of the ellipse defined by the points where the wavelet has zero 
value (zero-crossing), e = (1 — (c^/fTy) 4 ) 1 / 2 (for a y > a x ), and the sum 
s = a x + 0" 2 . It is alternatively described by the ratio of the semi-major and 
semi-minor axes of the Gaussian r = a x /a y , and the sum s = a x + a 2 . The 
axisymmetric Mexican hat is recovered for a x = <j v = 1, in which case r = 1 
(e = 0), and s = 2, and the normalization constant is unity, N(a x ,a y ) = 1. 

On the plane, the real Morlet wavelet is a typical example of a direc- 
tional wavelet. Its inverse stereographic projection on the sphere (see Fig. 
[3]) reads as (see also for similar projections): 



* (mor) M = \J^N{k) (l+tan 2 ^ 



cos 




with 7r -1 x = (2 tan(#/2) cos tp, 2 tan(#/2) sin tp) in Cartesian coordinates. 
The arbitrary wave- vector k = (k x ,k y ) controls the direction and the fre- 
quency of oscillation of the wavelet (k 2 = k 2 + k 2 ). The constant N(k) = 

(1 + 3e~ fc2 / 2 — 4e _3fe2 / 8 ) _1//2 stands for the normalization. Notice that for 
\k\ = 2, the real Morlet wavelet closely approximates at large scales to the 
second Gaussian derivative described in the following. 
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FIGURE 2: Mexican hat wavelet on the sphere for a dilation factor a = 0.4 
and different eccentricities. On the left, the axisymmetric Mexican hat: 
r = 1 (e = 0) and s = 2 (left). At the center and on the right respectively, 
the elliptical Mexican hat for r = 0.5 (e ~ 0.96825) and s = 2, and r = 
0.1 (e = 0.99995) and s = 2. Dark and light regions respectively identify 
negative and positive values. 




FIGURE 3: Real Morlet wavelet on the sphere for a dilation factor a = 0.4 
and a wave- vector k = (6,0) on the left, and for a dilation factor a = 0.4 and 
a wave-vector k = (2, 0) on the right. Dark and light regions respectively 
identify negative and positive values. 

The notion of filter steerability was first introduced on the plane \13\ 
[26] . and more recently defined on the sphere [29]. Just as on the plane, a 
directional filter ^ in L 2 (S 2 ,dQ) on the sphere is steerable if any rotation 
by X £ [0) 27r[ of the filter around itself R z {x)^ may be expressed as a linear 
combination of a finite number of basis filters Vl/ m : 



The weights k m (x)i with 1 < m < M, and M E N, are called interpolation 
functions. In particular cases, the basis filters may be specific rotations by 
angles Xm of the original filter: ^> m = R z (x m )^ '■ Steerable filters have a 
nonzero angular width in the azimuthal angle ip which makes them sensitive 
to a whole range of directions and enables them to satisfy the relation (|2.13j) . 
In the spherical harmonics space, this nonzero angular width corresponds to 
an azimuthal angular band limit N G N in the frequency index n associated 



M 




(2.13) 



m=l 
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with the azimuthal variable (p: 



for Inl > N. 



(2.14) 



Typically, the number M of interpolating functions is of the same order as 
the azimuthal band limit N. 

The derivatives of order in direction x of radial functions on the 
plane are steerable wavelets. The transfer of the steerability property (|2.13j) 
from the plane to the sphere is obvious since the inverse stereographic pro- 
jection is a radial operation, while the steerability only affects the azimuthal 
variable. The inverse stereographic projection of Gaussian derivatives there- 
fore give steerable wavelets on the sphere. They may be rotated in terms 
of M = + 1 basis filters, and are band- limited in ip at N = + 1. 
We give the explicit examples of the normalized first and second Gaussian 
derivatives. A first derivative has a band limit N = 2, and only contains the 
frequencies n = {±1}. It may be rotated in terms of two specific rotations 
at x = an d X = 7r /2, corresponding to the inverse projection of the first 
derivatives in directions x and y, ^> di and *ff a v respectively: 



R z (x) (w) = (u) cos x + (<*>) sin x- 



(2.15) 



The normalized first derivatives of a Gaussian (see Fig. 3} in directions 
x and y read: 



H 1 + tan 2 



1 + tan z 



tan — cos ipe 
2 Y 



tan — sin we 
2 ^ 



-2 tan 2 (0/2) 



-2 tan 2 (0/2) 



(2.16) 






FIGURE 4: First Gaussian derivative wavelet on the sphere for a dilation 
factor a = 0.4: from left to right, ^» a «(9 au ), \fr d y(9 au ) ; an( \ rotation by x = 7r /4 
oi ^d & {gau) Dark and light regions respectively identify negative and positive 
values [29]. 



A second derivative has a band limit N = 3, and contains the fre- 
quencies n = {0,±2}. It may be rotated in terms of three basis filters. It 
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reads indeed in terms of the inverse projection of the second derivatives in 



directions x and y, and respectively, and the cross derivative \& 



ds-d,-, 



as: 



R z ( X ) V°* (w) = (u) cos^ x+V 8 * M sin ^ X+^ ait>i M sin 2*. (2.17) 

The correctly normalized second derivatives of a Gaussian (see Fig. [5]) 
in directions x and y read: 



9| (gau) 



1PM 



1 + tan z - 1 - 4tan z - cos <£> e 



— ( 1 + tan z - ) ( 1 - 4tan z - sin z ip ) e 
4 



-2 tan 2 (0/2) 



-2 tan 2 (0/2) 



1 + tan — I I tan — sin 2ip ) e 



,-2tan 2 (6>/2) 



(2.18) 



z z z z 




FIGURE 5: Second Gaussian derivative wavelet on the sphere for a dilation 
factor a = 0.4: from left to right, ^ 9 l^) , $r 8 f(s a «) ; ^fyfoau^ and rota ti n 
by x = vr/4 of 1^1 (f a «). Dark and light regions respectively identify negative 
and positive values [29]. 



3. Directional correlation 

3.1 Directional and standard correlations 

The directional correlation (R^/\F) of a function F with a filter ^ is gener- 
ically defined as the scalar product of the signal with all SO (3) rotations of 
the filter [3'OJ. It therefore lives on 50(3), and explicitly reads in L 2 (SO (3), dp) 
as: 



(R(j>)V\F)= dn^* (R- l u))F{u) 



(3.1) 
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As discussed in subsection 12.11 if ^ is the specific dilation of a wavelet on 
the sphere, the directional correlation coincides with the wavelet coefficients 
of the signal, at the corresponding scale (see relation (|2.5p ). 

The standard correlation (Rq^\F) of F with is generically defined 
by the scalar product between the function F and the filter ^ translated 
at any point ojq = (#0,920) on the sphere, but for a fixed direction, i.e., a 
fixed value x = 0- The result of the standard correlation explicitly gives a 
square- integrable function in L 2 (S 2 ,dQ) on the sphere: 



The notation Rq simply denotes a three-dimensional rotation with x = 0. 
It distinguishes the standard correlation (Rq^\F) from the directional cor- 
relation (Rfy\F) when the arguments are not specified. 

Let us remark that, from relation (|2,13p . it explicitly appears that the 
directional correlation with a steerable filter ^ reduces to a M-terms linear 
combination of standard correlations with the corresponding basis filters 
\& m . In the particular case of an axisymmetric filter, there is no dependence 
at all of the correlation in the filter rotation \. The directional correlation 
with an axisymmetric filter is therefore strictly equivalent to the standard 
correlation. 

3.2 Pixelization and a priori computation cost 

The directional and standard correlations are defined for square-integrable 
functions on a continuous variable lo = (9, ip) on the sphere. The translations 
and rotations of the filter also form a continuous variable p = ({po,9o,x) S 
SO (3). Practical implementations are obviously based on a choice of dis- 
cretization for each of these variables, i.e., a pixelization of S 2 and £0(3). 
Let N p ~ (2-L) 2 represent the number of sampling points u in a given pix- 
elization of S 2 . The quantity 2L represents the mean number of sampling 
points in the position variables 9 and <p, or 9q and <po. A simple extrapola- 
tion of the Nyquist-Shannon theorem on the line intuitively associates L £ N 
with the band limit, or maximum frequency, accessible on that pixelization 
in the "Fourier" indices conjugate to 9 and (p for the signals and filters 
considered. For a sampling on 2L points in the direction x the same band 
limit L is associated with the conjugate Fourier index. Notice that in the 
wavelet formalism, the dilation parameter a G must also be discretized 
for practical purposes. 

Considering a simple quadrature, i.e., a discretization of the directional 
correlation integral, each scalar product on the sphere has an asymptotic 
complexity 0(L 2 ). The overall asymptotic complexity for the directional 
correlation (|3.ip . taking into account all discrete p = ((po,9o,x) on SO(3), 




(3.2) 
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is therefore of order C(L 5 ). We consider fine samplings of several megapixels 
on the sphere. To fix ideas, let us notice that the present NASA WMAP 
precision experiment on the CMB provides maps of the celestial sphere of 
around 3 megapixels. For a sampling associated with a band limit around 
L ~ 10 3 , the typical computation times for (2L) 2 multiplications and (2L) 2 
additions are of order of 1CP 2 seconds on a standard processor. We take this 
value as a fair estimation of the computation time required for a scalar prod- 
uct. Consequently, a unique 0(L 5 ) directional correlation would take several 
years at that band limit on a single standard computer. Moreover, depend- 
ing on the application, the directional correlation of multiple signals might 
be required. Typically, thousands of simulated signals are to be considered 
for a Monte Carlo statistical analysis. And for a wavelet analysis, multiple 
scales are to be considered for the filter. In conclusion, the directional corre- 
lation analysis of functions on the sphere is absolutely unaffordable for fine 
samplings with a band limit around L ~ 10 3 in 9, tp, and X- This conclusion 
remains when the use of multiple computers is envisaged. It is even strongly 
reinforced in the perspective of an analysis from finer pixelizations on the 
sphere. In particular, the forthcoming ESA Planck CMB experiment will 
provide 50 megapixels maps, i.e., L ~ 4 x 10 3 . 

The overall asymptotic complexity for the standard correlation, taking 
into account all discrete loq = (ipo, 6q) on S 2 , is of order C(L 4 ). On a single 
standard computer, the corresponding computation time through simple 
quadrature, at a band limit L ~ 10 3 , would be of the order of days. Such a 
calculation still remains hardly affordable, particularly when multiple signals 
and multiple scales are considered. 

3.3 Directional and standard correlations in harmonic space 

The Wigner D-functions coefficients {R^\F) mn of the directional correlation 
(R(p)^\F) living on SO (3) are given as the pointwise product of the spher- 
ical harmonics coefficients F[ m and ^J n . The following correlation relation 
holds: 

(R(p)y\F) = Y,^ E (mF) l mn Dtn(p), (3-3) 

leN \m\,\n\<l 

with 

(mF)mn = i^^nFlm- (3.4) 

Indeed, the orthonormality of scalar spherical harmonics implies the Plancherel 
relation (R^>\F) = XweN^| m |</ R^im-Fim- The action of the operator R(p) 
on a function G(u) in L 2 (S 2 , dtt) on the sphere reads in terms of its spheri- 
cal harmonics coefficients as: [R{p)G\ lm = ^|n|<z F) l mn {p)Gi n . Inserting this 
last relation for ^ in the former Plancherel relation finally gives the result. 



14 



Y. Wiaux, J. D. McEwen, and P. Viclva 



The standard correlation (R(ujo)^\F) lives on S 2 and could be decom- 
posed in its spherical harmonics coefficients. However, for non-axisymmetric 
filters, these coefficients do not appear as a simple pointwise product similar 
to (|3.4p . The easiest way to express the standard correlation in harmonic 
space is therefore to simply consider the relations (|3.3j) and (|3.4p with \ = 0. 

4. Fast algorithms 

4.1 Band-limitation 

The wavelet formalism defined in section [2] holds for any signal and any 
wavelet satisfying the admissibility condition (12, 8h . irrespectively of any 
band-limitation consideration. However, the band-limitation represents a 
necessary condition for obtaining precise numerical implementations. We 
therefore consider band-limited functions G at some band limit L E N on 
the sphere S 2 , i.e., G\ m = for I > L. From (|3.4p . the directional correla- 
tion of a band- limited signal F by a band-limited filter both with a band 
limit L on the sphere is thus also band-limited, with the same band limit: 

{m\F) l mn = for I > L. 

In practice, the signals F may generally be very precisely approximated 
as band-limited, through considerations relative to the physical data acqui- 
sition process. For the typical wavelets described in section [H ty a is also 
essentially band-limited, to very good approximation, provided that not too 
fine scales are considered (a 0). Under these conditions, the wavelet 
coefficients of (p, a) can therefore be calculated very precisely, or even 
exactly on equi-angular pixelizations, at suitable analysis scales. This is the 
scope of the fast directional correlational algorithms discussed in the next 
two subsections. 

We do not consider here the question of the signal reconstruction from 
its wavelet coefficients through formula (|2.7p . The corresponding numerical 
implementation would require an explicit discretization of both the scales a 
and the 50(3) variable p. First steps in that direction have been undertaken 
in [7]. 

4.2 Separation of variables 

The algorithm presented here for the directional correlation is based on the 
technique of separation of variables. 

The factorized form (12. ip of the spherical harmonics naturally enables 
one to compute a direct spherical harmonics transform by separation of the 
integrations on the variables 6 and (p. Conversely, an inverse transform may 
be computed as successive summations on the indices / and m, up to the 
band-limit L. Correctly ordering the corresponding operations provides a 
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calculation of direct and inverse spherical harmonics transforms in 0(L 3 ) 
operations [10]. This separation of variables for the spherical harmonics 
transforms may be performed on iso-latitude pixelizations on the sphere, i.e., 
pixelizations for which the sampling in 9 is independent of p, but where the 
sampling in <p may conversely depend on 9. This is the case for equi-angular 
pixelizations on the sphere. At a resolution L £ N, 2L x 2L equi-angular 
pixelizations are defined by a uniform discretization in 2L samples both for 
the angles 9 and p. Such a pixelization scheme defines pixels with areas 
varying drastically with the co-latitude [31] . In particular, on a 2L x 2L 
equi-angular grid, a sampling result on the sphere states that the spherical 
harmonics coefficients of a band-limited function with band-limit L may be 
computed exactly as a finite weighted sum, i.e., a quadrature, of the sampled 
values of that function [10] • HEALPix pixelizations (Hierarchical Equal Area 
iso-Latitude Pixelization) are also iso-latitude pixelizations, but where the 
sampling in ip explicitly depends on 9. Such a pixelization scheme defines 
l^Ysiefe pi xe l s °f exactly equal areas, for a resolution parameter N si( i e = 2 k 
with k G N. The computation of the spherical harmonics coefficients of 
a band-limited function is not theoretically exact on HEALPix grids, but 
can be made extremely precise by an iteration process [14] . These grids are 
notably used for the NASA WMAP CMB experiment and the ESA Planck 
CMB experiment. 

The very same reasoning based on the factorized form (|2.6p of the 
Wigner D-functions enables the calculation the inverse Wigner D-functions 
transform on SO (3) required by (|3,3p in 0(L 4 ) operations [T9j [20] . Con- 
sidering an iso-latitude pixelization for the angles 9 and ip on the sphere, 
the separation of variables for the Wigner D-functions transforms may be 
performed for any structure of the sampling in the third Euler angle x> 
potentially depending on 9 and (p. In particular, at a resolution parameter 
L £ N, one may consider a uniform discretization in 2L samples for x- Com- 
bined for example with an equi-angular pixelization at the same resolution 
for the angles 9 and (p on the sphere, this defines a 2L x 2L x 2L equi-angular 
sampling in p = (<p,G,x) on SO(3). 

Consequently, the algorithmic structure based on the separation of 
variables on iso-latitude pixelizations on the sphere may be summarized 
as follows. [18} 130] (a) Direct spherical harmonics transforms, and F/ m : 

0(L 3 ). (b) Correlation (R^\F) mn in harmonic space through d33|): 0(L 3 ). 
(c) Inverse Wigner L>-functions transform (R(p)^\F) on SO (3) through 
(|3.3p : 0(L 4 ). The global asymptotic complexity associated with the di- 
rectional correlation is thus reduced from 0(L 5 ) to 0(L 4 ) thanks to the 
separation of variables. For band-limited signals and filters, the numerical 
precision of the algorithm is simply driven by the precision of computation of 
the spherical harmonics coefficients. It is therefore very precise on HEALPix 
grids notably, and theoretically exact on equi-angular pixelizations. 
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4.3 Factorization of rotations 

The following algorithm for the directional correlation is based on the tech- 
nique of factorization of the three-dimensional rotations. 

The three-dimensional rotation operators R(p) on functions in L 2 (S 2 , dfl) 
on the sphere may be f actor ized as |251 [28l [22] 

R(vo,e ,x) =R(<Po- f R (o, - ,X + ~ ) • (4.1) 

The directional correlation relation (|3.3[) and the expression (|2.6p of the 
Wigner D-functions, matrix elements of the operators R(p), therefore give 
an alternative expression for the directional correlation of arbitrary signals 
F and filters on the sphere. We get indeed 

(R(p)9\F)= Yl (W> mm ^ i(mw+m " 0+nx) , (4-2) 

with the Fourier coefficients given by 

WU = e^- m W 2 £ d l m , m (|) d l m , n (|) % n F im , (4.3) 

i>c 

where C = max(|m|, \m'\, |n|), and with the symmetry relation d l , (9) = 

< •( <» m- 

For a band-limited signal F and a band-limited filter \& with band 
limit L £ N on the sphere one has |m|,|m'|,|n| < I < L. The factor- 
ized form of the imaginary exponentials enables the calculation of the in- 
verse three-dimensional imaginary exponentials transform required by (|4,2p 
in C(L 4 ) operations. Just as for the Wigner L>-functions transforms, con- 
sidering an iso-latitude pixelization for the angles 6 and (p on the sphere, 
the separation of variables for the three-dimensional imaginary exponentials 
may be performed for any structure of the sampling in the third Euler an- 
gle x- I n these terms, the directional correlation algorithm implemented 
on iso-latitude pixelizations for the angles and <p on the sphere through 
the factorization of rotations is structured as follows, (a) Direct spherical 
harmonics transforms, and F\ m : 0(L 3 ). (b) Correlation (R^\F) mm , n 
in harmonic space through f|4.3|) : 0(L 4 ). (c) Inverse transform (R(p)^\F) 
through (14. 2h : 0(L 4 ). The global asymptotic complexity associated with 
the directional correlation is thus also reduced from 0(L 5 ) to 0(L 4 ) thanks 
to the factorization of rotations E|. Again, for band-limited signals and filters, 
the numerical precision of the algorithm is simply driven by the precision of 
computation of the spherical harmonics coefficients. 



3 Notice that, while the Euler angles y?o and \ are i n the range tpo, x £ [0, the 
original range for 9q is 8q € [0, n], in order to cover the parameter space of SO (3). 
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4.4 Optimization with steerable and axisymmetric filters 

In terms of our rough estimations of subsection 13.21 the separation of vari- 
ables reduces the computation times on a standard computer from years to 
days for the directional correlation of band-limited signals and filters with 
band-limit L ~ 10 3 , typically sampled on megapixels maps. However, as 
already discussed, if a large number of simulations have to be analyzed, and 
at various scales of the filter, 0(L ) calculations remain hardly affordable 
even through the use of multiple computers. 

Steerable filters are typically considered with a small number of inter- 
polating functions M (see relation (j2.13fl ). that is also a small azimuthal 
band- limit N (see relation (|2.14p ) relative to L. The use of such steerable 
filters further reduces the asymptotic complexity for the directional corre- 
lation. On the one hand, the directional correlation with a steerable filter 
^ reduces to a M-terms linear combination of standard correlations with 
the corresponding basis filters *$> m . For M <C L, the asymptotic complexity 
of a directional correlation reduces to that of a standard correlation, with 
an a priori 0(L 4 ) complexity, to which is simply added the 0(L 3 ) linear 
combination which arises from (I2.13|) . On the other hand, on iso-latitude 
pixelizations on the sphere, either the technique of separation of variables, 
or the factorization of three-dimensional rotations can be applied to the 
standard correlation, by setting x = in the relations (|3.3p or (|4.2p respec- 
tively. For a steerable filter with a small azimuthal band limit N <C L, the 
Fourier index n, with \n\ < N, can be excluded from asymptotic complexity 
counts. It readily appears that the corresponding asymptotic complexity 
for the two algorithms hence reduces to 0(L 3 ), on iso-latitude pixelizations 
on the sphereQ At L ~ 10 3 , our rough estimation of computation times is 
reduced from years to tens of seconds. This renders the computation easily 
affordable, even when multiple signals and multiple scales are considered. 

Details on the algorithmic structure, computation times, memory re- 
quirements, and numerical stability of the corresponding implementations 
on HEALPix and equi-angular grids on the sphere may be found in [22] for 
the factorization of rotations, and in [30J for the technique of separation of 
variables and the optimization with steerable filters. Notice in that regard 



Considering also 9o £ [0, 2it[ puts the result on the parameter space of the three- 
torus T 3 , which covers twice the parameter space of SO(3). In that context, the 
relation (|4.2|l is understood as a three-dimensional inverse Fourier transform, which 
can be calculated in 0(L 3 log 2 L) operations on a 2L x 2L x 2L equi-angular grid on 
50(3) by the use of the standard Cooley-Tukey fast Fourier transform algorithm. 
This optimization however does not reduce the overall asymptotic complexity for 
the directional correlation, still driven by the computation of (|4.3|) in 0(L 4 ) oper- 
ations. 

Let us remark that the issue of the sampling in \ is not relevant for steerable 
filters. The proper rotations by \ £ [0, 2n[ are indeed analytically driven, and thus 
with infinite precision, by the relation (|2.13[) . 
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that a further optimization of the algorithm based on the separation of vari- 
ables and with steerable filters may be achieved on equi-angular pixelizations 
on the sphere. It relies on the fact that Wigner L>-functions transforms may 
be decomposed into linear combinations of spherical harmonics transforms, 
which therefore drive the overall asymptotic complexity for the directional 
correlation. On 2L x 2L equi-angular pixelizations, these spherical harmon- 
ics transforms may be computed in 0(L 2 log 2 L) operations through the 
Driscoll and Healy algorithm [lOl HI US], if the associated Legendre poly- 
nomials are pre-calculated. As discussed above, the sampling theorem on 
equi-angular pixelizations on the sphere also renders the calculation exact. 

The axisymmetry of a filter A{6) on the sphere is an extreme case of the 
steer ability, for an azimuthal band limit N = 1: A[ n = for |n| > 1. In that 
case, we already emphasized that the proper rotation by x nas no effect on 
the filter and the directional correlation reduces to a standard correlation. 
At each scale, the wavelet coefficients of a signal with an axisymmetric filter 
therefore live on the sphere S 2 rather than on 5*0(3). The directional cor- 
relation relation (|3.4p consequently reduces to the following standard form, 

giving the spherical harmonics coefficients {RoA\F) lm of the correlation of 
a signal F with an axisymmetric filter A as the pointwise product between 
the filter's Legendre coefficients Ai, and the spherical harmonics coefficients 
of the signal Fi m : 

(ihA\F) lm = 2-KA\F lm . (4.4) 

The correlation of a band-limited signal with a band-limited axisymmetric 
filter (A[ = for I > L) is therefore readily computed in the harmonic 
space of S 2 . On iso-latitude pixelizations on the sphere, the direct spherical 
harmonics transform of the signal, and the inverse spherical harmonics trans- 
form of the correlation, can simply be computed by separation of variables 
in the spherical harmonics. This provides an algorithmic structure with 
0(L 3 ) asymptotic complexity, which again can be reduced to Oil? log 2 L) 
on equi-angular pixelizations. 

5. Conclusion 

A new field of complex data processing has emerged in many areas of science. 
Scalar and tensor data, often distributed on nontrivial manifolds, come up 
at continually increasing resolutions. Powerful signal analysis techniques 
need to be developed to process such datasets. 

In this paper, we first reviewed recent formal developments for the con- 
tinuous wavelet decomposition of signals on the sphere. Second, we detailed 
advances in the definition of the corresponding fast directional correlation 
algorithms. 

These generic developments can find many applications in various fields 
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such as computer vision (omnidirectional cameras, ...), biomedical imag- 
ing (functional magnetic resonance imaging, ...), geophysics (signals on the 
Earth's surface, ...), or astrophysics and cosmology (signals on the celestial 
sphere, ...). In that regard, the important results already obtained in cos- 
mology through the wavelet analysis of the cosmic microwave background 
strongly illustrate the fact that the formalism developed represents a pow- 
erful tool for complex data processing on the sphere [23J . 
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